function omega = cgns_vorticity(ubase,grid,udofs,uuu,varargin)
% omega = cgns_vorticity(ubase,grid,udofs,uuu)
% omega = cgns_vorticity(ubase,grid,udofs,uuu,disc)
%
% Project the vorticity field on the velocity base. To simplify
% the projection, we use a simple algebraic averaging.
%
% If the input parameter disc is present and different from 0, then
% the reconstruction is discontinuous, otherwise a continuous
% reconstruction is returned.

 % precompute the local projection matrix
 for i=1:ubase.pk
   for j=1:ubase.pk
     M(i,j) = sum( ubase.wg .* ubase.p(i,:) .* ubase.p(j,:) );
   end
 end
 Mi = inv(M);
 % local projector
 for i=1:ubase.pk
   wgp(i,:) = ubase.wg .* ubase.p(i,:);
 end
 uMiwgp = Mi * wgp;
 clear M Mi wgp

 if(nargin>4)
   disc = varargin{1};
 else
   disc = 0;
 end
 if(disc)
   omega = zeros(ubase.pk,grid.ne);
 else
   omega = zeros(udofs.ndofs,1);
   nnn = omega;
 end
 for ie=1:grid.ne
   % vorticity at the Gauss nodes
   for id=1:grid.d
     idx = grid.d*(udofs.dofs(:,ie)-1) + id;
     uue = uuu(idx);
     for jd=1:grid.d
       gue(jd,:,id) = uue' * permute(ubase.gradp(jd,:,:),[2,3,1]);
     end
     % coordinate transformation
     gue(:,:,id) = grid.e(ie).bi' * gue(:,:,id);
   end

   gomega = (gue(1,:,2) - gue(2,:,1))';

   if(disc)
     omega(:,ie) = uMiwgp*gomega;
   else
     idx = udofs.dofs(:,ie);
     omega(idx) = omega(idx) + uMiwgp*gomega;
     nnn(idx) = nnn(idx) + 1;
   end
 end

 if(not(disc))
   omega = omega./nnn;
 endif

return

